Stroboscopic wave packet description of 
time- dependent currents through ring-shaped 

nanostructures 

Martin Konopka^ Peter Bokes* 

"^Department of Physics 
Institute of Nuclear and Physical Engineering 
Faculty of Electrical Engineering and Information Technology 
Slovak University of Technology in Bratislava 
Ilkovicova 3, 812 19 Bratislava, Slovakia 

^European Theoretical Spectroscopical Facility 
(ETSF, www.etsf.eu) 

Received: date / Revised version: date 

Abstract 

We present an implementation of a new method for explicit simula- 
tions of time-dependent electric currents through nanojunctions. The 
method is based on unitary propagation of stroboscopic wave packet 
states and is designed to treat open systems with fluctuating num- 
ber of electrons while preserving full quantum coherence throughout 
the whole infinite system. We demonstrate the performance of the 
method on a model system consisting of a ring-shaped nanoj unction 
with two semi- infinite tight-binding leads. Time-dependent electron 
current responses to abrupt bias turn-on or gate potential switching 
are computed for several ring configurations and ring-leads coupling 
parameters. The found current-carrying stationary states agree well 
with the predictions of the Landauer formula. As examples of gen- 
uinely time-dependent process we explore the presence of circulating 
currents in the rings in transient regimes and the effect of a time- 
dependent gate potential. 

PACS numbers: 73.63.Rt 
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1 Introduction 



Transition from the state with zero electric current to a state with nonzero 
current in nanoelectronic devices is a technologically important example of 
non-equilibrium quantum dynamics of an open many-electron system. It ex- 
hibits several significantly different time scales related to relaxation of elec- 
trons with important consequences for the prospects of their applications 
in nanoelectronics. The common problem in theoretical quantum transport 
description in both stationary and time-dependent situations is proper in- 
clusion of the semi-infinite leads. The number of degrees of freedom in the 
nanojunction, which are explicitly considered, has to be sufficiently small in 
order to be numerically tractable. Hence the boundary conditions at the 
junction must be correctly described and included. An approach widely used 
for stationary systems are the non-equilibrium Green's functions (NEGF) 
with lead self-energies [U |2] which include the physics induced by the leads. 

The description of boundary conditions of the finite nanojunction be- 
comes more complicated in time- dependent problems. The non-stationary 
quantum transport has been addressed for non- interacting electrons [3] , within 
the time-dependent density-functional theory (TDDFT) [H [5] or even the 
many-body perturbation theory [6]. It has been shown that self-energies - 
now time-dependent - could in principle be used also in this lone as 

the leads are treated at the TDDFT level. In practice, it is a quite com- 
plicated approach especially for realistically structured systems so that sim- 
plifications of the formalism are necessary [7J |S]. Chen et al. proved the 
so-called holographic electron density theorem and developed a new method 
for time-dependent open electronic systems [9]. Further effort in this direc- 
tion led to computationally more efficient density-functional tight-binding 
method [S]. Very recent works of Chen et al. on time-dependent quantum 
transport are based on Liouville-von-Neumann equation for single-electron 
density matrix [10J. Another recent work by a different group [H] uses gen- 
eralised master equation approach to mesoscopic time-dependent transport. 
Non-equilibrium thermodynamical theory of interacting tunnelling transport 
has been presented by Hyldgaard (12] . 

Another class of methods to tackle time-dependent transport with open 
boundary conditions has been in development [131 ttH US] • The scope of these 
methods is wider than just elastic transport. The methods are know as cor- 
related electron-ion dynamics (CEID). The dynamics is based on Ehrenfest 
molecular dynamics and its extensions. The electronic degrees of freedom in 
CEID are divided between the central system and its environment which is fa- 
cilitated by the formalism of one-particle density matrices. A damping term 
has been introduced used in the open-boundary equations of motion (13] ; 
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this term keeps the environment close to a reference state. In further de- 
velopment of the method [14], open boundaries have been introduced in a 
new way which represents an explicit realisation of an external battery. We 
remark that this method does not conserve coherence between injection and 
subsequent scattering of the electrons. 

In the present work we address the open-boundaries time-dependent quan- 
tum transport problem using the stroboscopic wave packet (SWP) basis 
method, principles of which has been provided in the works JTBJ [TTj . We study 
time-dependent currents through model-system nanojunctions formed by 
small rings inserted to electric circuit formed by mono-atomic leads. We also 
present new developments of the stroboscopic wave packet approach (SWPA) 
which are necessary in order to describe the systems with atomistic structure. 
While in the orginal works [TBI [TTj structureless electrodes and very simple 
tunnelling barriers have been used to demonstrate the method, here we work 
with atomistic models of the electrodes and more complex nanojunctions. 
The formulation presented here solves the time-dependent Schrodinger equa- 
tion (SchE) for independent electrons within the tight-binding (TB) approx- 
imation. SchE is solved numerically with state vectors expanded in the stro- 
boscopic wave packet basis representation (SWB) [THl [TT]. The SWPA has 
been designed to be specifically suited for time-dependent transport through 
nanojunctions. Its main advantage is to make possible explicit integration 
of equations of motions of nanoj unction degrees of freedom with full quan- 
tum coherence preserved throughout the whole infinite system. This can be 
viewed also as an open-system treatment with correct inclusion of boundary 
conditions. In the present paper we focus our study on simple model cases 
of rings and monoatomic leads with possible extension to realistic systems. 
Our main interest are transient currents developed in response to abruptly 
applied bias and time-dependent gate potential. We also refer to stationary 
results obtained by other authors using analytical methods. 

Transport properties of atomic-scale sized rings have been studied by 
several authors in most cases on the Hiickel/tight-binding level of descrip- 
tion and in the stationary regime. Effect of the asymmetric position of lead 
on transmittance has been studied in Ref. [18] by means of Green's func- 
tions. Authors of work [Tj5] have developed a source-sink potential method 
for convenient analytical treatment of molecular electronic devices including 
conjugated systems. This approach has been used in Ref. [20] to obtain the 
form of the transmission with explicitly given dependence on the molecular 
skeleton and its connection to leads. Authors of Ref. [21] generalised the 
so-called waveguide approach [22] (WGA) to systems described by the TB 
approximation. Resulting methodology - the tight-binding waveguide ap- 
proach (TBWGA) - can be interpreted as a route to an exact solution of the 
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stationary SchE for given class of TB models. An explicit formula for trans- 
mittance has been derived [21] for rings composed of identical atoms and 
identical nearest-neighbour couplings and coupled to two leads with equal 
chemical potentials in both of them. The authors discussed conditions for 
occurrence of circulating currents in such quantum rings. More recently 
Sparks et al. formulated the stationary problem more generally, considering 
a multibranch device in a TB approximation [23J . They provided an explicit 
formula for system transmittance valid for a large class of systems including 
those with different chemical potentials in different leads. The formula has 
been obtained by an exact solution of the stationary SchE for the whole sys- 
tem and has also been related to Green's function analysis. Explicit results 
in Ref. [23] facilitate the study of quantum interference effects in stationary 
regime for wide class of systems. We refer to these results especially when the 
long-time behaviour of our time-dependent method is discussed. An inspiring 
computational work was done by Saha et al. |24j. The authors study a true 
atomistic model of a quantum-interference-controlled molecular transistor 
formed by an 18-annulene attached between zigzag graphene nanoribbons. 
The device was studied by a multi-terminal NEGF-DFT formalism [25J and 
the current-switching effect was confirmed. 

On the experimental side, ring-shaped tunnel nanoj unctions can be formed 
by cyclic organic molecules attached between two electrodes (see for exam- 
ple [26j[271|28l[29]). Conductors can be formed by metallic (usually Au or Pt) 
electrodes or by graphene nanoribbons [301 EI] ■ These works employ direct 
attachment of particular molecules in between metallic electrodes (Pt-C or 
Au-C bonds) using mechanically controllable break junctions. Such contacts 
are reported to be more conductive and stable than previously more com- 
mon junctions employing a bridging thiol group (thus forming a metal-S-C 
bond) [32]. Another realisation of nanometer-scale sized rings is fabrica- 
tion of structures on proper substrates [331 EH ES] . This can be done using 
molecular-beam epitaxy, wet etching and optical or electron-beam lithog- 
raphy. Finally we note that the concept of a quantum interference driven 
transistor started to be more explicitly discussed in literature relatively re- 
cently [36]. It still represents an experimental challenge and a major bunch 
of experimental results is yet expected to come. 

The paper is organised as follows. In Sec. [2] we describe the model of the 
atomic ring. In Sec. [3] we explain the implementation of the stroboscopic 
wavepacket method. In Sec. [4] we provide the formula which we implement 
for electron current calculations. Sec. [5] describes stationary results as a basis 
from which we move into non-stationary regime in section [6j which contains 
our main results. 



4 



2 The model of the ring with contacts 



We consider a linear chain of atoms, each two being a lattice constant a 
apart, with one finite ring which presents an obstacle for the flow of electrons 
(Fig. [T]). All couplings are considered within the TB approximation in which 
we limit our treatment to one orbital per atom The chain is periodic apart 
from the ring region. Hence, the TB Hamiltonian is of the form 

H(t) = H° + H 1 + H 2 (t) (1) 

where 

oo oo 

H° = 2J e a \ a l + ^2 ^b(4+i q / + a \ a i+i) ■ ( 2 ) 

l=— 00 l=— 00 

a\ and ai are the fermionic creation and destruction operators, e is a constant 
on-site energy in the chain, and £3 is the TB hopping parameter here assumed 
to be negative. The important simplifying assumption is that there is only 
one localised atomic orbital (state \l)) per each site. The H° represents the 
unperturbed chain, later referred to as the lead's Hamiltonian. 

H 1 is a stationary, periodicity-breaking term. Its form corresponds to 
a TB ring structure, indicated in Fig. [I] (Specified values of N and n are 
taken only as an example there.) In the simplest case, the modifications to 
the periodic chain Hamiltonian H° due to the ring presence will be expressed 
by the following matrix elements of H l : 

H-N+1,N = ~~ ^ B 

H N,1 = *B • (3) 

H-N+l,n = 

All other matrix elements of H l are zeros. 

Applied bias will be modelled by time-dependent shifts of the on-site 
energies. If a bias voltage U(t) is applied then the on-site energies of the left 
lead will take values e + eU(t), e being magnitude of electron charge. The 
operator H 2 (t) expressed in the atomic-orbital basis has the following matrix 
elements: 

{ell(t) 5i )V , I < (in the left lead) 
\eU{t) 5 l)V , 1 < I < N (within the ring). (4) 
, all other I, V 

The on-site energies within the ring, if not differently specified, will thus 
take values e + eU(t)/2. The on-site energies in the right lead will remain 
unchanged, equal to the equilibrium value e. 
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Figure 1: Topology of the studied systems and the numbering of atoms. 
Throughout the paper we use iV for the total number of the ring atoms as 
well as for the index of the atom just above the left vertex. We use n to 
denote the index of the right vertex atom. Symbols M and M' are provided 
for convenience to display relation to the notation of Ref . [21] ■ 



Three exceptions from the ab ove pr escription of the perturbation V(t) 



H 1 +H 2 (t) are studied: (i) In Sec 



6.1.2 



where a uniform slope of the potential 



energy within the ring is used instead of the spatially constant value eU(t)/2, 

(ii) in Sec. 6.3, where a time- and branch- dependent gate potential is is used, 

(iii) finally consideration of variable coupling strength between the ring and 
the leads in Sec. 16.41 



3 The stroboscopic basis set 



The stroboscopic wavepacket basis set consists of wavepackets (Fig. [2]) con- 
structed from the eigenstates of the unperturbed (lead's) Hamiltonian 
The mathematical expression for the basis set vectors is [TBI Ej 



\n, a, m; t) = exp 



--(mr n + t)H° 
n 



^£ n 



d£ \8, a) 



(5) 



where H° is the lead's Hamiltonian and \£,a) are its eigenstates normalised 
so that (S,a\S', a') = 5(8 — £') 5 aja >. Each basis function or wavepacket ^ 
is uniquely characterised by three indexes: the band index n, the time shift 
index m, the degeneracy index a and time t. 
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The band index n. The unperturbed Hamiltonian ^ gives the disper- 
sion relation 

E(k) = e + 2t B cos(ka) . (6) 

The corresponding energy range [e + 2t B ,e — 2£ B ] is in the SWPA divided 
into non-overlapping bands [15]. n th band has its energies £ G [£ n _i,£ n ]. In 
our present work we use two equally wide bands [37] spanning the whole TB 
energy range. The band index n then takes values 1 a 2. 

The time-shift index m. Within each band, different, mutually or- 
thogonal basis functions are obtained by time shifts mr n , where the time 
step 

T - = as: (7) 

is set by the energy width of the band A£„ = £ n — £ n ~i- rn attains all integer 
values, but in numerical calculations it is restricted torn = — m max , . . . , +m max 
as it is discussed at the end of this section. 

The degeneracy index a. Apart from the energy, each eigenstate of 
the lead's Hamiltonian has further quantum numbers. In the present system, 
the only one is the direction of propagation of the Bloch eigenstates: those 
propagating from the left to the right (index a = +1) and the opposite ones 
(index a = —1). 

The time t. The time t in the notation indicates that the basis state 
is unitarily propagated [TTJ by the lead's Hamiltonian (J2]). The use of this 
"moving" basis set significantly simplifies the form of the SchE in the SWB 



representations (see Eq. 11 below), very much in the spirit of using the inter- 
action representation in formal perturbation theory in quantum theory. At 
each given time t the SWB vectors form an orthonormal system since 

(n, a, m; t\n', a', rri; t) = 5 n>n >5 a!a >8 mim > . (8) 

The set of basis states ^ is complete if m max is infinite. 

Given the form of Hamiltonian 0, each electron in the system evolves 
independently as described by the SchE 

ih^Mt)) = R(W{t)) ■ (9) 

State vectors \^(t)) for each electron are expanded in this basis set: 

2 ^max 

n=l a=±l m=— m max 
4(2m max +l) 

^ Mt)\o;t), (10) 

0=1 
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propagation direction 








lattice site index 



Figure 2: Schematic view on stroboscopic wave packets \n, a, m;t) at time 
t = 0. Actual picture shows absolute values of their projections to atomic 
orbitals |Z), i.e. the quantities \(l\n, a, m; t) | as functions of the lattice site 
I. Band index is chosen to be n = 1, propagation direction a = 1 and 
the m indices run through the range —8, —7, . . . , +8. The very low value of 
m max = 8 is chosen for convenient visualisation and is actual only for this 
scheme. 



where we have introduced a composite index o = (n,a,m). The SchE pro- 
vides a linear system of differential equations for each of the involved elec- 
trons. Due to the time-evolution of the basis functions given by the operator 
H°, the equations for amplitudes A (t) include matrix elements of the per- 
turbation only: 

(\A (t) 4 ( 2mmax+1 ) 
ih -jr= E (o;t\V(t)\o';t) A ,(t) , (11) 

o'=l 

with V(t) = H 1 + H 2 {t). This is of key importance for our method: basis set 
wavepackets, which are localised far from the central region feel either zero 
perturbation (in the right lead, and hence their matrix elements are zeros) 
or they feel only the spatially uniform bias-induced potential of the left lead; 
see eq. (|4]). For such matrix elements we have 

(o;t\V(t)\o';t)neU(t)S 0>0 , (12) 
in the left lead, for packets far from the central region 

meaning that probability amplitudes of the wavepackets evolve freely, (inde- 
pendently of other amplitudes) far in the left lead. Situation is even simpler 
far in the right lead where V(t) is zero and corresponding amplitudes do not 
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evolve at all. The constant or freely evolving amplitudes need not be explic- 
itly included into simulation and this omission in principle does not have any 
impact on accuracy of the simulation. In practice, since the SWPs are not 
strictly localised, by using the finite cutoff m max we introduce certain error 
into simulations. 

The simple form of equations of motion (11 ) has been possible due to the 
employment of the moving (i.e. unitarily propagating) basis set. The sim- 
plicity is in that the matrix elements are computed only from the interaction 
term V(t) = H 1 + H 2 (t) of the total Hamiltonian Q. 

Wavepackets corresponding to vectors ^ unitarily propagate in time 
which implies that each particular packet will travel away from the central 
region after some time. The result would be that (for finite m max ) at large 
times the region would not be covered by basis set wave packets at all and 
no electrons would be presents in the central region at long times. In detail, 
the unitary propagation of the wavepackets is such that during a time inter- 
val r n [defined below ([5])] each packet moves exactly to the position of the 
neighbouring packet (either on its left or on its right, depending on the prop- 
agation direction). If we had a very large basis set (m max — > oo) and were 
looking at a movie of the wavepackets only at the "stroboscopic" times m'r n , 
we could not see any motion of the packets. For any finite m max we could 
observe similar static picture only for a limited time and in a limited spatial 
region because the wavepackets constantly propagate from one's position to 
the other's. 

To prevent the gradual disappearance of the basis functions from the 
central region in simulations, we periodically insert new basis wave packets 
into the system at every period r n . The newly inserted packets are localised 
far away from the centre, but moving toward it. To avoid the increase of 
the total number of basis vectors, we also periodically remove basis wave 
packets (those which has left the central region). In other words, we remove 
the leading wave packet from the train of the propagating packets and at- 
tach a new wave packet at the tail of the train. (See also Fig. |2}) This 
removal/insertion procedure is accomplished for each band (indexed by n) 
independently; different bands could in principle have different widths and 
consequently different parameters r n . For given band n there are two such 
procedures: one for a = — 1 (packets propagating from right to left), the 
other for a = +1 (packets propagating from left to right). In this way we 
keep the basis set vectors in the region of interest (the "central region" ) and 
also keep the number of the vectors constant. 

It is necessary to have m max sufficiently large so that the basis vectors 
|n, a, m; t) with m « m max are decoupled from other basis vectors (see more 
reasoning below) and hence the removal of the basis vector |n, a, m max ; t) 
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does not affect the quantum dynamics of the electron, amplitude of which 
has been removed. 

In addition to the basis vector removal/insertion procedure, we also track 
individual electrons in the simulation. Each newly inserted basis state, if 
belonging to normally occupied band (n = 1 in this work), is set to be occu- 
pied at the instance of the insertion. Also, an electron, which through time 
evolution gets far from the central region, is removed from the simulation ac- 
cording to a proper criterion. Therefore the simulation explicitly describes an 
open system in which the total number of electrons generally fluctuates [38J . 

The initial conditions applied to the probability amplitudes are such that 



■^■n,a,m (0) 



1, for n 
0, for n 



1 

2 



for all a, m, 



(13) 



i.e. only the lower band is initially occupied. 

Matrix elements of V{t) are evaluated using the overlaps (l\o;t). In Ap- 
pendix we derive that the overlaps can be expressed by the formula 



(l\n, a, m; t) = (l\o; t) 
i 




e (mr n + 1) 



aKA — 2% (mT n + t) cos /C 

h 



d/C 



(14) 



Quantities K. n = k n a are dimensionless wavenumbers corresponding to ener- 
gies 8 n through the TB dispersion relation (|6]). The cutoff / max on atomic 
sites indexed by I has to be chosen sufficiently large in order to cover all 
SWPs ^ included into simulation. For majority of calculations we use 
m max = 352 with rather conservatively chosen Z max = 3016. In one case we 
use m max = 1200 and corresponding Z max = 7576. Integrals in (14) have to 
be evaluated numerically. 

The SWPA as described above and used in this work does not include any 
electron-electron (e-e) interactions except for keeping them always in strictly 
orthogonal states. We expect that e-e interactions with effects like Coulomb 
blockade among others become more important for weakly coupled rings 
(although they can never be neglected in any accurate quantitative descrip- 
tion) [39]. Further developments of the SWPA will include also realistic e-e 
interactions at least on a mean-field level of description like in time-dependent 
Hartree-Fock theory or in TDDFT. Presently we could only include a simple 
dynamical mean-field interaction according to a chosen model. Having done 
this we have not found any significant interesting impacts of the interaction. 
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Hence we decided not to include any such results in the present work as it 
would present an unnecessary complication. 

The system of the equations (11) is solved using the modified- midpoint 
method |40j. Its accuracy and stability is fully sufficient for given problem. 
As a complementary treatment applicable for stationary currents, we com- 
pute exact stationary currents for given TB Hamiltonian ([I]) in cases when 
it is time-independent. See Sec. [5] for more details. 



4 Electron current formula 



In this section we briefly provide an expression for local electron current. 
For this purpose we consider one-dimensional infinite TB chain, specified in 
Sec. |2j For such a system the local current operator at bond between sites 
lo and Iq + 1 is given by the formula (see Ref. p], p. 162 therein). 



l s ,l 



1 

ih 



Hi +i,i alj Q+1 a s j - Hi j Q+1 a\^a Sj i 0+1 



(15) 



with s — ± 



1 



being the spin index. Taking the expectation value in a single-electron state 
\^{t)) we obtain 



I$L(t) = ^H l0+hl0 (t) (¥(f)|/ + 1) (l |*(f)> + cc. . 



ih 



(16) 



Superscript (1) marks that it is a current caused by single electron; we must 
sum up over all electrons in the system to obtain the total current. If we 
express the state vector in the stroboscopic basis orbitals then we obtain 



ih 



J2A* o (t)(o;t\l + l) 



o=l 



,o'=l 



+ CC. 



(17) 



The factor of 2 has been added to take into account two electrons differing 
only by their spins. As for the sign convention used for the current expressed 



by eqs. (16), (17), it indicates flow of particles (electrons), not the charge. 



The same purely local result for I; (i) as shown above could be derived 



also for currents in particular ring arms. Hence we use formula (17) to 



compute time-dependent current (contribution from one electron) through 



11 



any chain (lead or arm of the ring) of the complete system. The total current 
is obtained by summing up contributions (17) generated by all individual 
explicitly included electrons in the system: 



(18) 



ele=l 



5 Stationary currents through ring nanojunc- 
tions 

Currents through TB rings have been studied in literature mostly in sta- 
tionary regimes. It is convenient to compare quasi- stationary currents from 
our method to exact stationary results for the TB model. For this purpose, 
the whole system - ring with leads - is assumed to be composed of identical 
atoms described by the simple TB Hamiltonian with all couplings equal to tn. 
The only perturbation to periodic chain Hamiltonian H° are then terms (pi) 
which represent a perturbation to the topology of the linear chain. It was 
shown in Ref. [2T] that in the stationary regime there are both conducting 
and insulating configurations depending on where the leads are attached to 
the ring. Specifically, all odd-numbered rings are always conductive irrespec- 
tive of the choice of attachment site. On the other hand, even-numbered rings 
can exhibit both behaviours. These findings are confirmed and extended by 
the exact TB results of Ref. [23] • 

We have obtained stationary results from exact eigenstates of the TB 
Hamiltonian H° + H 1 + H 2 defined in Sec. [2] This approach is, for the Hamil- 
tonian used in our work, equivalent to the method described in Ref. [23J. As 
specified in Sec. |2j the applied bias U is modelled by lifting on-site energies 
in the left (source) lead at the value e L = e + eU. The on-site energies in the 
right lead (drain) are kept at e H = e. On-site energies of all ring atoms are 
kept at the intermediate value e s = e + eU/2 = (e L + e R )/2. TB couplings 
are unchanged, i.e. they all are equal to £b- The eigenstates for such system 
are expressed in the TB basis, 



They depend on the eigenenergy E and the direction of propagation a. For 
brevity we do not show these dependences explicitly. We compactly represent 



00 




(19) 



/=— 00 
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the amplitudes by composite formula 



with 



V /(/C S Z) + Q /(-/CgZ) 

.A /(/C L Z) + B /(-/C L Z) C /(/C R Z) 

Tf(ic s i) + gf(-)c s i) 



/(zcz) 



P{E) m 
2tt 



(20) 



(21) 



The specific form of (20) is an abbreviated non-standard notation mimicking 



the spatial positions of particular system chains and has to be understood 
in the sense that ipi = A f{K.\J) + B /(— JC^l) for Z 6 left lead, ipi = 
T /(/CsZ) + £ /(— /CsZ) for / G lower ring branch, etc. We introduced the 
dimensionless k-number 

K = ka (22) 

with a being the lattice constant. p(E) is the local density of states (DOS) 
of an ideal infinite chain for a particular branch of the system and for the 
TB model is expressed as 



p{E) 



(23) 



e, identified with e L , e s or e R , is the value of the on-site energies in given 
branch of the system. In the model under study it also represents Fermi 
energies of particular bulk systems. Indices L, S and R stand for the left 
lead, small system (the ring) and the right lead, respectively |41j. When 
convenient we use also notations pl, Pr and ps to distinguish DOSs in the 
left lead, right lead and in the small system (ring). Parameters A and B 
describe wavefunction of the left lead, parameter C in the right lead, J 7 and 
Q in the lower ring arm and finally V and Q in the upper ring arm. 

An electron being in particular eigenstate is then interpreted as par- 
tially reflected and partially transmitted, with transmittance computed as 



T(E) 



(24) 



This simple formula is applicable also in cases when the Fermi levels in the 
two leads are different. In such a case we need additional factors involving 
the respective group velocities. These factors have been explicitly inserted 



into functions (21) hence are not explicitly present in the formula (24). The 
results from exact eigenstates for the ring of iV = 16 sites and n = 7 are 
shown on Fig. [3j 
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Figure 3: Transmittances of the ring depicted in Fig. [I] (N = 16, n = 7 or 
equivalently n — 11) calculated from exact eigenstates of the TB Hamilto- 
nian; see text and Ref. [23]. The plots with higher e L values are vertically 
shifted for convenience and the shifts are visually enhanced by horizontal 
dashed lines. Parameters e L are on-site energies in the left lead. Right leads 
has always set e R = 0. The small system (ring) has e s = (e L + e R )/2. The 
transmittances are computed and shown only for energies at which extended 
propagating eigenstates exist. £b < is the tight-binding hopping parameter, 
magnitude of which is used as the unit of energy throughout this work. 



2e 
~h 



Having transmittances available, we compute the stationary electron cur- 
rents in the leads using the Landauer formula 

/ T(E)dE . (25) 
'e f 

We remind that the transmittance T(E) depends parametrically also on the 
chemical potentials in the leads (which are equal to e L and e H ) and on the 
on-site energies of the ring atoms which are all equal to e s . In the eigenstate 
approach we always use e R = and e L = ell (applied bias) and the on-site 
energies of the ring are set to e s = (e L + e R )/2. 



The integral in Eq. (25) is done numerically. In Fig. [4] we show the 
comparison of the exact stationary currents and the currents obtained by 
the time-dependent wavepacket approach. Results from the SWPA (discrete 
symbols) represent long-time quasi-stationary values from time-dependent 
calculations. We see overall agreement between the two approaches, es- 
pecially for low biases. The results from eigenstates are exact for given 
Hamiltonian. The SWPA [TBI E] is i n principle exact too, but its accu- 
racy is lowered by the finite cutoff to the number of basis functions which is 
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— exact eigenstates 




U/(\tJ/e) 



Figure 4: Stationary electron currents vs. bias through the ring system 
depicted on Fig. fl] (N — 16, n — 7) calculated by the two different ap- 
proaches. Solid black line shows electron current obtained from the Lan- 
dauer formula p5| in which the transmittances, shown on Fig. [3j have been 
computed from the exact solution of the TB model. Discrete symbols show 
long-time quasi- stationary values from our time-dependent SWPA at several 
basis sizes given by the values of m max [HJ [TF] . In most of this work we use 
either m max = 352 or m max = 1200. In these and all other numerical results 
throughout this work we consider the periodic tight-binding model of the 
leads with zero Fermi energy in equilibrium state. Out of equilibrium, Fermi 
level in the left lead is lifted up to the bias value U while keeping the right 



lead at zero Fermi energy. See main text, most importantly subsection |6.1.2 
for further details. 



2iVbands(2m max + 1). The impact of the finite cutoff is typically not important 
at low biases but increases with the bias and the obtained accuracy is then 
lower at high bias conditions. Convergence of the current with the cutoff 
becomes very slow at large m max . Several examples computed with different 
basis set sizes are shown on Fig. |4j The steady-state current in the leads 
is always underestimated [12] when finite m max is used. To understand the 
finite basis set impact, we first remark that the dynamics of an individual 
electron becomes (quite obviously) more accurately described when using 
larger m max , hence providing more accurate current (17) generated by single 
electron. Second, the increased m max results in a larger number of electrons 
iV e i e included into the total current formula (18). As a consequence, the pre- 
cise dependence of the current on m max is a complicated function. The much 
slower convergence at higher values of m max can qualitatively be understood 
as a result of larger spread of the SWPs at hight m (shown on Fig. [2]). In 
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the limit of m — > oo the packets become fully delocalised. Individual high-m 
packets contribute very weakly to the total local current at given point in 
space. The less converged results in absolute terms are pronounced espe- 
cially at larger biases. The convergence is worse also in relative terms: at 
U = 0.2 \t B \/e we reach 95.5% of the exact limit (Fig. g|. At U = 0.5 we get 
91.6% and at U — 1.8 only 89.3%, all at m max = 352. The worse convergence 
at larger biases can be qualitatively understood on the basis of the equations 
of motion (11) and bias-related matrix elements (12). Using a finite m max 
in (11) we drop some portion of the exact equations of motion, in particular 
some of the terms (12) which describe the interaction with bias-induced lifts 
of on-site energies in the left lead. Using higher bias increases significance 
of those matrix elements and hence makes convergence more difficult. The 
finite basis set size would not impact the accuracy of our approach if the 
stroboscopic basis states with large \m\ indices were well localised in space. 
In reality, even wavepackets with very large |m| indices exhibit long tails 
towards the central point of the lattice where the nanojunction is placed. 
Ideally, having basis wavepackets (those with large |m| indices) with no tails 
in the nanojunction region would result in perfectly converged results. (Such 
an ideal situation is impossible because of the dispersion which spreads the 
wavepackets. Our further effort in development of the SWPA is directed to 
overcome the convergence issues despite the presence of the dispersion.) 



6 Simulations 

6.1 Bias voltage switch on ring-shaped nanojunction 

Time-dependent bias is modelled by given time-dependent (but spatially uni- 
form) lift of the TB on-site energies in the left lead only. The atoms in the 
right lead have their on-site energies all set to zero, which is also the value 
of Fermi energy in equilibrium. The atoms within the ring have their on- 
site energies set to the half of the bias value applied at given time although 
we test also different model in subsection 16.1.21 In our simulations we first 
"equilibrate" given system by letting it evolve at zero bias (and zero temper- 
ature) for time up to t sw = 500 7i/|£b|- The "equilibration" evolves the non- 



interacting many-electron initial state of the system, described by eq. (13), 
into (again non-interacting) stationary many-electron state of the Hamilto- 
nian with the localised perturbation (the ring system and eventual variations 
of some of the on-site energies or interatomic couplings). We have verified 
that the equilibration period is sufficient to obtain converged transient elec- 
tric currents in the sense that these currents are practically independent of 
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Figure 5: Upper panel: Local electron current as a function of the two 
variables, lattice site and time, computed for conductive ring configuration 
with N = 1Q and n = 7. Lattice sites from 1, . . . , 16 belong to the ring 
structure. The applied bias of 0.5 |te |/e was abruptly switched on at time 
t sw = 500 | . Currents through vertex sites (here 1 and 7) do not have 
well defined values and certain interpolated values have been used for visual- 
isation purposes. Lower panel: Analogous plot but now for insulating con- 
figuration with n = 8. Basis set size for these plots uses m max = 352. Time 
is expressed in units of ^/|£b| and local electron currents (on the z axes) in 
units of e\t-Q\/h. On this figure we do not apply any artificial smoothening 
of the data (see caption to Fig. [6] for example). 

particular choice of t sw if t sw is at least 500 ^/ |te | - The transient electron 
currents are evaluated according to formula ( fi~7| ), with applied summation 
over all contributing electrons in the system. We study how these transient 
currents depend on several ring parameters like the ring size, placement of 
the terminals and the magnitude of applied bias. Typical behaviour of the 
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transient electron currents is shown in Fig. [5j It is noticeable that both the 
temporal and spatial evolution must be considered. The intersect of the two 
ridges corresponds to the time when the bias was turned on and to the spatial 
location of the ring. The slope of the ridge corresponds to the Fermi velocity, 
t>F = (l/h) dkE(k = tt /a) = —2at-B/h, so that the time interval between the 
instance of switching-on the bias and the arrival of the density or current 
perturbation to a site I is given by the expression t\ = la/vp. The ridge on 
the left side corresponds to the reflected electron impulse propagating back 
into the left lead. Similarly, the other ridge represents the transmitted wave- 
front in the right lead. The function values are mostly positive meaning that 
the currents in both leads are positive (i.e. electrons flow from the source 
lead to the drain lead). Results for several specific cases are discussed in the 
following subsections. Presented ring sizes N vary from 16 to 18. We have 
calculated also several other sizes, starting from N = 5, not shown here. 
Plotted local currents will typically correspond to a lattice site positioned 
far from the ring, in most cases site / = 120, which is displaced about 100 
lattice constants from the considered rings. 

6.1.1 Dependence on the location of drain vertex 

The dependence of the time-dependent current on the position of the right 
(drain) lead is shown in Fig. [6j We show how and if the transient currents 
(induced by the abrupt bias switch at time t sw = 500 ^/|£b|) depend on cho- 
sen drain vertex site n. In Fig. [6] we plot transient local currents obtained 
from our simulations for three subsequent rings sizes and several drain-lead 
attachments. These rings represent various different kinds of electronic trans- 
port properties which can be found in rings. The local currents shown have 
been calculated at lattice site 120. The finite basis set adds rapid oscillations 
on the computed local current^ We smooth out the oscillations by averag- 
ing them over the time interval of about r. This averaging does not modify 

1 The period of these artificial oscillations is always r, the quantity defined by eq. d7b, 
here equal to iTh/\ts\ (same for both the bands). The oscillations are larger at sites closer 
to the ring. All this is explained by the localised nature of the stroboscopic basis states 
(the wavepackets) and the finite number of them. The train of the wavepackets (see Fig. [2]) 
travels in real space in such a way that during a time interval [t, t + t] each packet places 
itself exactly to the position which was occupied by its neighbour at time t. Since we use 
the finite cutoff on the basis set, the quality of the description at given point of space 
fluctuates in time with the period of r. Instead of using a huge value of 77i max (which 
would be prohibitively expensive) we smooth out the oscillations by averaging them over 
the time interval of about r. Typical timescales of the processes studied in the present 
work are in most cases much larger than r and the relevant effects are not affected by the 
oscillations or by their artificial smoothening. 
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Figure 6: Transient currents computed in the drain lead at site I = 120 for 
several different drain lead attachments (vertices specified by the numbers 
n = 5, 6, 7, 8 shown at individual plots) of the rings composed of iV = 16, 17 
and 18 atoms. In all cases the applied bias U = 0.5 I^b | /e has been switched 
on abruptly at time t sw = 500 ^./|t B | ■ Basis set size again uses m max = 352. 
The plots have been smoothed by taking running averages over an interval of 
about 7ih/\t-B\. The dashed lines represent analytically computed stationary 
currents as obtained by the exact-eigenstate method described in Sec. [5j 

essential features of the plots including the characteristic times and rates. 
Example of the raw unaveraged data are shown in Fig. [7| The "no slope" 
black curves in that figure are obtained in the same system as plots n = 7 
and n = 8 in Fig. [6j apart from the time averaging. 

The onset of the electronic responses appears at times around 550 VI*b|, 
which is delayed about 50 time units after t sw , caused by the distance from the 
ring to the observation site / = 120 [ti ~ (I — N)a/vp ~ 50 fo/\tB\ }■ We see, 
that similarly as in the stationary regime, the time-dependent simulations 
also exhibit the two very different regimes - the conducting and insulating 
one. However there is the peak of transient current even in the insulating 
configurations after the bias is switched on. This feature is shown in the 
lower panel of Fig. [5] and will be discussed also below. 

N — 16 rings. The left panel of Fig. [6] shows results for rings of size 
N = 16, but differing by the position of the drain vertex sites which run 
through the sequence n = 5, ... ,8. The odd-numbered cases {n = 5 and 
n = 7) correspond to conducting configurations. The well conducting state of 
the ring is given by constructive interference (CI) of the electron amplitudes 
in the two ring branches. CI are most significant in the case of symmetrically 



19 



attached ring N — 16, n — 9, the configuration with equally long branches 
(not shown in the figure). As we can see from the figure, and in agreement 
with former theoretical analysis [2"T| 123], significant CI is possible for several 
vertex configurations. The two local currents (n = 5 and n = 7) reach simi- 
lar long-time limits and the time needed to build up the current is the same 
for each drain vertex site: At ~ lOft/liel- The dashed lines in Fig. [^repre- 
sent stationary currents as obtained from the exact-eigenstate approach. At 
intermediate times, around 600 time units, the dynamically computed cur- 
rents are in good agreement with the stationary approach. The agreement 
would become slightly worse at larger times when finite basis set errors take 
effect. See also Fig. [4] for comparisons of stationary currents also at other 
biases. On the other hand there are the even-numbered cases which show 
almost zero currents after a short transient effect. The peaks of the tran- 
sient currents (the plots with n = 6 and with n = 8) are very similar each 
other. The blocking status is again due to quantum interference, now the 
destructive one. Results for other values of n, not shown in plots, also con- 
firm that the transient characteristics are practically the same for all drain 
vertex attachments within the particular group (conductive or insulating). 

N — 17 rings. Central panel of Fig. [6] shows results for rings with size 
N = 17. In this case there is no such a distinct separation into conductive 
and insulating configurations, the N = 17 rings are all conductive. Neither 
constructive nor destructive inteferences are now perfect. However, similarly 
as in the insulating configurations and contrary to the conducting configu- 
rations of even-sized rings, we now observe peaks in transient currents. In 
addition there is a substantially longer relaxation, now lasting for about 100 
time units. The peak current depends on the drain vertex site in an oscil- 
latory manner. (More DC-conductive rings in this class exhibit lower peak 
currents.) N = 17 rings have shown to be more difficult from computa- 
tional point of view. The effect of the incomplete basis set would show up 
at long-time stationary currents which would be, by estimate, 10-25% lower 
compared to exact results from eigenstates. 

N = 18 rings. In this case (right panel of Fig. [6J we again obtain two 
distinctly different behaviours. The difference from the N = 16 ring is that 
now the odd-numbered drain vertices correspond to insulating configurations. 
The long-time limit of the time-dependent treatment would relate to exact 
stationary results similarly as for N = 17 rings. In addition, the ring with 18 
atoms and with the right vertex at site n = 5 exhibits pronounced circulating 
currents as will be discussed in separate subsection |6.2 

The results show, that for even-numbered rings, the transient effect is 
very short, lasting for about 10 time units. Surprisingly, the odd-numbered 
ring iV = 17 exhibits longer relaxation, lasting for about 100 time units, 
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Figure 7: Transient local currents through the system with N — 16 ring 
computed in the drain lead (again at site I = 120) for one conducting (n = 7) 
and one insulating (n = 8) configuration and for two models of on-site energy 
setting within the ring: the model with spatially constant profile ( "no slope" , 
black line) and the model with linearly varying energies ("slope", green or 
grey line). The abruptly turned-on bias has the magnitude 0.5 |te | /e. Basis 
set size for these plots are the same as on Fig. [6j In contrast to other figures, 
the shown currents were not averaged over short time scales. 

although the initial development upon the bias switch is equally rapid as for 
N = 16 and N = 18 rings. The long-time asymptotic values of currents from 
our simulations exhibit either conducting or insulating character which is in 
many cases in a good agreement with the analytical stationary results. 

6.1.2 Role of potential-energy profile within the ring 

The time evolution of the system in our model is considered in the independent- 
electron approximation. The on-site energies within the ring are all set to 
the half of the instant bias value. The on-site energy profile within the ring 
is spatially constant. One may expect that the results would be different if 
we had considered spatial variation of the on-site energies within the ring. 
The energies might, for example, be computed on-the fly as dependent on 
the electron charge density which would be a simple dynamical mean-field 
model. Another model for ring on-site energies may be their linear variation 
(slope) from one terminal of the ring to the opposite terminal such that the 
on-site energies of the vertex atoms reach those in the corresponding leads. 
We checked both these models and found that they do not have any signifi- 
cant impact on the electron dynamics in the tested cases. This is illustrated 
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in Fig. [7] for a ring composed of 16 atoms and the drain vertex being ei- 
ther at site n = 7 (conducting configuration) or at site n = 8 (insulating 
configuration). As we can see the differences between the spatially uniform 
and the spatially varying model are essentially negligible for both conduct- 
ing and insulating configurations. This finding may not be universally valid. 
However, to keep the models in the present work simple, we stick at the 
spatially uniform profile of on-site energies within the studied rings, with 



obvious exceptions of an applied gate potential (see Sec. 6.3). 



6.2 Circulating currents 

It has been discussed in several studies (e.g. [21J) that circulating currents 
(CCs) may arise in certain ring configurations. These analysis were done 
for stationary currents. In this subsection we study circulating currents in 
transient regime upon the abrupt switch of the bias. By definition, a CC 
at given instant of time in a ring structure (like that on Fig. [T]) exists when 
the currents in the two ring branches have opposite directions in the sense 
that, for example, the electrons in the shorter branch flow from the left 
terminal to the right one while the electrons in the longer branch flow from 
the right terminal to the left one. The currents in the ring branches are 



calculated using formulae (18) and (17). These formulae represent strictly 



local currents: h (t) is in fact the current between sites Iq and lo + 1, i.e. a 
current through the bond. Our calculation of the current in given ring branch 
in addition involves averaging over the interatomic bonds of given branch. 

Inspection of the currents calculated from exact transmittances of Ref. [23] 
shows that stationary CCs are possible for many ring sizes. Actual occur- 
rence of the CCs depends also on chosen drain site and on applied bias. Of 
the studied rings we found most pronounced CC in N = 18 ring with its 
drain vertex at site n = 3 and n = 5. The results for the n = 5 structure are 
shown in Fig. [8] 

Upper panel in Fig. [8] shows the exact stationary electron currents, cal- 
culated with the approach described in [23] (lines), alongside with the quasi- 
stationary values from our simulations computed with the basis-set parame- 
ter m max = 1200 (isolated symbols) [13]. Red (dark grey solid line) and green 
(light grey dashed line) plots display currents in the two ring branches and 
are relevant in verifying if CC is present. The sign convention used on the 
figure is such that currents in both branches have positive signs for the flow 
of the particles from left to right, i.e. state without CCs. The negative value 
of one of the currents indicates the presence of the CC in the ring. From the 
upper panel of Fig. [8] we see that such a negative electron current flows in 
stationary regime in the longer ring branch for a wide range of voltages (0 
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Figure 8: Electron currents in the system lead and in the two ring branches, 
computed for N = 18, n = 5 ring. Upper panel: Line plots show stationary 
currents computed from exact transmittances [23] • Point symbols show quasi- 
stationary values obtained in our simulations in a long-time limit (times 
about 1000 units). Basis set corresponding to m max = 1200 (larger than 
in other figures) has been used. Lower panel: Time-dependent electron 
currents induced by abrupt bias switch at time t sw = 500^/^1, computed 
for bias 0.694 |^b I / e - This value of bias corresponds to the largest circulating 
current in the stationary regime. Colour and line-style coding is the same as 
on the upper panel. Currents in the branches are spatial averages over sites 
within particular branch (excluding the vertex sites). Currents in the right 
lead (computed at site 120) have been smoothened similarly as described in 
caption to Fig. [6] m max = 1200 has been employed also here. 



to 0.78 |tB | /e) . Most pronounced CC is found at a bias of U = 0.694 |te|/e. 

Lower panel of Fig. [8] extends the results at the bias 0.694 |£b|/c into the 
non-stationary regime. The non-stationarity is due to the abrupt bias switch 
at time t sw = 500 h/\t^\. The colour and sign convention is the same as for 
the upper panel. Negative currents show up in the longer ring branch. Black 
plot on the lower panel of Fig. M represents current computed at site 120 (i.e. 
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the usual current in the right lead at a position 115 sites beyond the right 
vertex atom)|^| 

Finally, the dip in the right-lead current (black plot on the lower panel of 
Fig. [I]) at times ranging from about 630 to 730 h/\t-Q \ is caused by the numer- 
ical cutoff on the stroboscopic basis set size. This unphysical dip tends to 
be more pronounced in certain ring configuration, for certain quantities and 
for high biases. Ring configuration (N,n) = (18,5) is significantly affected 
by the finite basis set error at certain time interval. This is also the reason 
why we used parameter m max = 1200 to obtain results shown in Fig. [8] while 
for most of other figures we use m max = 352. The increased cutoff improves 
accuracy of time-dependences especially at intermediate times. 

6.3 Time-dependent currents controlled by gate field 

Quantum interference effects between electrons in the two ring arms can be 
controlled by an applied electric field. The field value can be time-dependent 
in general. In our approach it is modelled by a spatially uniform lift of on-site 
energies in one of the branches of the ring (not including the two terminal 
atoms); the affected atoms have indices n + 1, n + 2, . . . , N in the indexing 
scheme of Fig. [T} 

6.3.1 Abrupt gate turn-on 

In our first studied application of the gate potential, we control the electron 
transport using an instantaneously switched gate field which is applied to 
the two N = 16 rings differing by their drain terminal indices n. The field is 
turned on at time £ g = 800 h/\t^\, after the system has been evolving under 
the constant bias of U = 0.5 {tsl/e so that it is essentially in a stationary- 
current regime before the time when the gate field is lifted. The simulations 
give results which are depicted in Fig. [9j The local current is evaluated at 
the lattice site 120. We observe that the gate field induces transient effects 
lasting for about 10 — 40 time units. The transient effects depend rather 
strongly on particular bias values. The application of the higher gate field is 
interesting; it would allow the ring to work as a field driven switching device, 

2 The smoothcning of the artificial oscillations (see previous footnote) is not always 
perfect as can be seen from the bottom panel of Fig. [HJ especially the black and red plots. 
This is a purely technical issue caused by a too large interval between recorded data from 
the simulation and by the width of the window used to compute the running averages. It 
is interesting that the green plot in the graph does not show the oscillations. In fact they 
are present but are very small in magnitude. The impact of the finite basis set on the 
presence of the oscillations is often quantitatively different in different parts (leads and 
ring branches) of the whole system. 
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Figure 9: Local electron currents computed at site 120 when gate field is 
abruptly turned on at time t g = 800/*i/|£b|- The applied bias is constants 
in these plots and equal to U — 0.5 |^b | / e. Values of the gate potential 
are shown in legends and use the same units as the bias. The ring has 
N = 16 atoms. The plots in the left panel use drain vertex attached at site 
n = 7 (conducting configuration) while the plots in the right panel use n = 8 
(insulating configuration). The dotted lines represent analytically computed 
stationary currents as obtained by the exact-eigenstate method described in 
Sec.[5j The exact stationary currents are even functions of the gate potential. 

i.e. a nanoscale-sized transistor, changing the conducting ring into non- 
conducting or vice versa. The switching is a consequence of effective change 
in Fermi wavelength of electrons in the controlled (longer) branch and hence 
changing the interference pattern which shows up in the resulting current. 
The dynamically computed currents at longer times agree well with the exact 
values from eigenstates although differences arise especially at higher gate 
potentials. As in other calculations, the source of the differences is the finite 
cutoff on the number of stroboscopic basis packets used. 

6.3.2 Harmonic gate potential 

Another option to control the current flow is to employ a harmonically os- 
cillating gate potential. For this purpose we decided to use larger rings with 
sizes N ranging from 70 to 90. Such rings provide longer characteristic times 
of temporal changes. We study the following three (N, n) configurations: 
(70,61), (80,71), and (90,81). All the three rings thus have one longer and 
one shorter branch. The gate field is again applied to atoms with indices 
in the range n + 1, . . . , N, i.e. the atoms of the shorter branch now. The 
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Figure 10: Main graph: Average lead currents as functions of the sinusoidal 
gate potential for three different ring structures: iV = 70, n — 61 (black line 
with circles), N = 80, n = 71 (red line with squares), N = 90, n — 81 (green 
line with diamonds). The bias is always U = 0.5 |£b|/c (time-independent). 
The gate field acts on the atoms n + 1, . . . , N. The amplitude of the gate 
potential is always V g o = 0.25 Inset: Detailed time-dependences 

shown for three particular points from the red plot of the main graph: the 
case of N = 80 with uj g = 0.08 (dark blue plot, solid line), 0.12 (brown plot, 
solid line), and 0.16 |£b|M (orange plot, dashed line). Technical and basis set 
parameters used for these results are the same as for most of the other plots; 
in particular, m max = 352. 



constant bias is again U = 0.5 |ts | /e. The magnitude of the sinusoidal gate 
potential is always V g o = 0.25\t^\/e. The potential is harmonic with various 



angular frequencies as shown on Fig. 10 The potential starts to act at time 



t g = 800^/ltel- The inset of Fig. 10 shows typical time dependencies of the 
current in leads as obtained for the case (N, n) = (80, 71). The three plots of 
the inset have been recorded for three different gate angular frequencies u g . 
0.08 (dark blue, solid line), 0.12 (brown, solid line) and 0.16 of |£b|/^ (orange, 
dashed line). The corresponding periods T g = 2tt /uj g are approximately 78.5, 
52.4, and 39.3 of h/\t-B,\- The current response is generally quasi-periodic but 
anharmonic. (The initial transient effect after the gate is turned on is not 
discussed here.) The temporal dependence of the current is given by an in- 
terplay of the harmonic gate potential and internal ring effects given also 
by its size. At low gate frequencies u g (for example 0.08 |te|/^ or less, see 
the blue plot, dark solid line, in the inset), the current oscillations typically 
(not always) exhibit periodic quasi-harmonic pattern at the twice of the gate 
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frequency. The doubled frequency arises from the symmetric dependence of 
the system transmittance on the gate potential. (V g and — V g have the same 
effect on stationary currents as could be seen from the analytic formulae of 
Ref. |23j or calculated by formalism of Sec. |5j) Rapid driving is not followed 
by the current in this sense as can be seen from the inset. For example, the 
gate field with co g = 0.16 |£b|/^ (the orange plot with dashed line on the fig- 
ure) or at higher gate frequencies results in current oscillations at the same 
frequency. Intermediate driving frequencies (e.g. 0.12 |£b|/^> the brown plot, 
solid line, on the figure) yield periodic but anharmonic time dependences. 

It is interesting to compute time averages from the oscillating lead cur- 
rents. The three plots in the main graph of Fig. 10 show that the average 
current for given ring structure depends on the gate frequency (while the 
bias and the gate amplitude are kept unchanged). The maxima and minima 
of the average currents vary with varying ring structure: larger rings have 
the extrema shifted towards lower frequencies. Inspection of the plots show 
that these variations fulfil the law 5u/u = —5N/N which is expected from 
elementary considerations about resonant frequencies. In this way we have 
an evidence that the oscillatory pattern of the average current plotted in 
Fig. 10 arises from the internal ring resonances. Given the nanometer-scale 
size of such ring structures, the characteristic resonant frequencies lie in the 
optical domain and we do not investigate higher gate frequencies that those 
shown in the graph. 

Electric current dynamics could be investigated also for a fixed ring size 
N and varying drain terminal n. However, our inspection has shown that 
the dynamics is more interesting (i.e. the average currents exhibit more pro- 
nounced oscillations as functions of u g ) when the source and drain terminal 
are relatively close to each other. 



6.4 Reduced coupling to the leads 

Through previous sections it was assumed that all nearest neighbour cou- 
plings were identical along the whole composed system, i.e. including the 
small system (the ring). The couplings between the leads and the ring were 
then quantified by the hopping parameter £b < (same as in the lead and 
in the ring). In this section we study a more general situation when the 
couplings within the leads again remain equal to ts and the same parameter 
is used also to bound the atoms included in the ring. The couplings between 
the leads and the ring will however be characterised by a new parameter 
^bs < (lead-system coupling). 

We will report on regimes with |£bs| — |^b|- In the limiting case of 
I^bsI — > the leads and the system (the ring) become decoupled and the 
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Figure 11: Local electron currents evaluated at site 120 for various lead- 
system coupling strengths £bs- Legends indicate values of tss relative to 
the magnitude of the lead's hopping parameter t-g. As in previously shown 
results, bias U = 0.5 |te |/e is abruptly turned on at time t sw = 500/VI^b|- 
Number of ring atoms is iV = 16 in both panels and drain is attached to 
sites n = 7 (left panel) and n = 8 (right panel). The dotted lines in the left 
panel (conducting configuration) show stationary currents computed from 
the exact-eigenstate approach. 



internal energy levels of the ring attain discrete values. In intermediate cases 
of < |£bs| < |^b| we expect broadening of the ring levels. The internal ring 
electronic structure will have an impact on the transport characteristics. 



An illustration of our findings is shown on Fig. 11 The time-dependent 



currents have been computed for the rings of N = 16 atoms, left panel 
corresponding to drain vertex site n = 7 (conducting configurations) while 
right panel shows plot for n = 8 (insulating configurations). In all cases the 
bias U = 0.5 |£b|/c was abruptly switched on at time t sw = 500 /x/ |^b | • The 
magnitude of the transient effect as well as the currents (in the conductive 
system) behave as expected regarding to the lead-system coupling strength. 
The insulating ring retains its property also at reduced couplings. However, 
at intermediate coupling range tes £ [—0.50, —0.25] there are decaying os- 
cillations in the electron current persisting for a significant period of time, 
thus extending the transient effect to a period of about 100fVI^B|- Inspec- 
tion shows that angular frequency uj of these decaying oscillations is related 
to the applied bias by the relation hu = elf/2, in agreement with previous 
analytical [T7] as well as numerical p] results. 
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7 Conclusions 



Recently proposed stroboscopic wave packets [TBI E] have been developed 
to be applicable to systems employing explicit atomistic level of modelling. 
Time-dependent transport of electrons in an open system with localised per- 
turbation was studied. The localised perturbation had a ring structure de- 
scribed in the tight-binding approximation. Such a system is of interest due 
to quantum interference effects that affect its transport properties. It might 
serve as a field driven quantum interference current switching device [36] 
- a nanoscale-sized transistor. We have demonstrated the potential of our 
newly developed method based on the unitarily propagating stroboscopic 
wave packets to describe open systems with fluctuating number of explicitly 
included electrons while the whole system has an infinite number of elec- 
trons which can be neglected in a good approximation. In our method full 
quantum coherence is preserved throughout the whole infinite system. The 
method can be used for systems with one- dimensional semi-infinite leads. 
It is capable of spatially resolved description of transient effects like those 
caused by an abrupt bias switch of a gate field application. However, its 
most useful application would be found in systems which are exposed to 
long-lasting varying external fields or biases and for which it is important 
to preserve quantum coherence in the description. For short-time simula- 
tions one could use a less expensive model with cyclic boundary conditions. 
However, such an approach would become prohibitively expensive for large 
simulation times as it would require to consider many explicit atoms to de- 
scribe the leads. In contrast, our method does not have any such limit on 
simulation time. The weakness of the present implementation of the stro- 
boscopic basis set description is slow convergence of relevant results with 
increasing number of basis functions. The finite basis set demonstrates itself 
in a transient unphysical drop of lead current in some of the studied systems. 
Depending on the studied configuration, long-time stationary currents may 
also sometime be significantly underestimated. Other consequence of the fi- 
nite stroboscopic basis set are small artificial rapid oscillations in computed 
quantities which can however be smoothed out and usually does not prevent 
us from capturing relevant physical time-dependent effects. In our ongoing 
work we will consider a generalisation of the basis set in order to reduce the 
finite-basis set errors and to reduce the number of basis states needed. 
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A Derivation of the (Z|n, a, m;t) overlaps 



Here we derive formula (14) for the overlaps between the TB atomic orbitals 
\l) and the SWPs \o;t) = \n,a,m;t). In the SWPA [HI H7j, the eigenstates 
\£ , a) of the basis-generating Hamiltonian H° [given by eq. (p2J) in this work] 



use the energy "normalisation" condition which is 

(£, a\£', a') = S(£ - £')8 a>a , . (26) 
See also sec. [3] and eq. ^ therein. This condition leads to the form 

\£, a) = —== life, a) (27) 

\dS\ 
V \dk\ 

where \k, a) is the usual Bloch wave in one dimension with the k-number of 
magnitude k and the propagation direction labelled by a = ±1. Normalisa- 
tion of the Bloch waves is assumed to be 

(k, a\k', a') = 5{k — k')5 a ^ a > . (28) 

(The k-numbers will be restricted to the interval k = [0, ir/a] where a is the 
lattice constant.) In case of the TB model with the set of the orthonormal 
atomic orbitals |/) (one orbital per atom) we obtain formula 



1 / 1 

\£,a) = —== J— V e iaKl \l) (29) 

' - /2|* B |sin/C V 2vr U y } 



with the summation running over all lattice sites. /C = ka is the dimen- 
sionless wavenumber. The eigenstates \£, a) are used to construct the SWPs 
according to eq. Now we use that formula and write down the overlaps 
in the form 



(l\n, a, m: t) = — r / exp 



' (l\£,a)d£ (30) 



i 

■- (mT n + 1)£ 
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where we have utilised the equation 



H°\S,a) =S\£,a) 



(31) 



Using the expression (29) for the eigenstates and substituting it into eq. (30) 
leads to an integral expression with integration variable S. With the aid of the 
TB dispersion relation (JS1), the integral can be transformed to the integration 
variable /C. The final form of the expression for the overlap (l\n,a,m;t) is 
then provided by formula (14). We evaluate these integrals numerically. 
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